Astronomy & Astrophysics manuscript no. 13666 


© ESO 2009 


December 25, 2009 





Letter to the Editor 



o 
o 

o 

Q 

w 
6 

o3 



> 

oo 

00 

On 
O 



3D simulations of supernova remnants evolution 
including non-linear particle acceleration 

Gilles Ferrand 1 , Anne Decourchelle 1 , Jean Ballet 1 
Romain Teyssier 1 ' 2 , Federico Fraschetti 3 4 

1 Laboratoire AIM (CEA/Irfu, CNRS/INSU, Universite Paris VII), CEA Saclay, bat. 709, F-91191 Gif sur Yvette, France 

2 Institute for Theoretical Physics, University of Zurich, CH-8057 Zurich, Switzerland 

3 LUTh, Observatoire de Paris, CNRS-UMR8102 et Universite Paris VII, F-92195 Meudon Cedex, France 

4 Lunar and Planetary Lab & Dept. of Physics, University of Arizona, Tucson, AZ, 85721, USA 



draft 2.8 - Dec. 22, 2009 



ABSTRACT 



If a sizeable fraction of the energy of supernova remnant shocks is channeled into energetic particles (commonly identified with 
Galactic cosmic rays), then the morphological evolution of the remnants must be distinctly modified. Evidence of such modifica- 
tions has been recently obtained with the Chandra and XMM-Newton X-ray satellites. To investigate these effects, we coupled a 
semi-analytical kinetic model of shock acceleration with a 3D hydrodynamic code (by means of an effective adiabatic index). This 
enables us to study the time-dependent compression of the region between the forward and reverse shocks due to the back reaction of 
accelerated particles, concomitantly with the development of the Rayleigh-Taylor hydrodynamic instability at the contact discontinu- 
ity. Density profiles depend critically on the injection level n of particles: for r\ < 10~ 4 modifications are weak and progressive, for 
77 ~ 10" 3 modifications are strong and immediate. Nevertheless, the extension of the Rayleigh-Taylor unstable region does not depend 
on the injection rate. A first comparison of our simulations with observations of Tycho's remnant strengthens the case for efficient 
acceleration of protons at the forward shock. 
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1. Introduction 

Supernova remnants (SNRs) are believed to be the major con- 
tributors to Galactic cosmic rays. But although acceleration of 
electrons in these objects is well established, direct evidence of 
the acceleration of protons remains elusive: recent detections of 
non-thermal gamma rays do not yet allow unambiguousl y dis- 
tingu ishing leptonic and hadronic contributions (see e.g. Gabici 
2008 and references therein). 

To assess particle acceleration in SNRs, a promising alter- 
native approach consists in diagnosing the impact of energetic 
particles on the SNR dynamics. Indeed, if SNRs are efficient 
accelerators, as claimed, then energetic nuclei must make a siz- 
able impact on their e volution (see e.g. iJones & Ellison] 1 1991. 
Malkov & Drury 2001). Such a study is now possible thanks to 
the performance of modern X-ray observatories such as Chandra 
and XMM-Newton, which allow the s tructure of young Ga lactic 
SNRs to be probed in great detail. IWarren et al.l (120051) have 
been able to determine the positions of the three waves (forward 
shock, contact discontinuity, reverse shock) in Tycho's SNR with 
rather good accuracy, and have found that it does not match any 
pure hydr odynamic model. The same effect has been repo rted in 
SN 1006 (Cassa m-Chenai et al.ll2008l iMiceli et alJl2009h and in 
Cas A (Patnaude & Fesen 2009). This is evidence that not all the 
kinetic energy from the explosion is converted into heat down- 
stream of the shock, but that a sizeable part is channeled else- 
where - probably into energetic particles. 

This effect has been studied by iDecourchelle et all ([2000), 
using ID self-similar simulations c oupled with a simple mode l 
of non-linear acceleration (from Berezhko & Ellison 1999). 



They have shown how the shocked region shrinks in the case of 
efficient acceleration of particles at the shocks: as the injection 
fraction rises the shocks get closer to the contact discontinuity. 
These results have been confirmed and extended by ID hydrody- 
namic simulations of radially symmetric SNRs (see lEllison et al.l 
2007 and references therein). 

But the contact discontinuity between the shocked ISM and 
the shocked ejecta is known to be hydrodynamically unstable: 
this interface is subject to the Rayleigh-Taylor instability, as the 
ejecta are being decelerated by the ambient medium of lower 
density. Thus the morphologi cal study of a SNR requires a 3D 
(or at least 2D) modelling (see |Chevalier et al ]|1992llDwarkadasl 



2000, IWang & Chevalier! 1200 ll and references therein). As the 
instability develops the contact discontinuity is profoundly mod- 
ified: the ejecta protrude inside the shocked ISM, forming char- 
acteristic finger-like structures. These features are particularly 
apparent in Tycho's SN R, where th e ejecta exhibi t a fleecy aspect 
in bot h X-rays (Wa rren et al.l2 005) and in radio (Velazque z et al.l 
1998). Thus, to dia gnose the back reaction of energetic particles, 
one needs to take h ydrodynamic instabilities into account. To as- 
sess their impact, iBlondin & Ellison] (1200 ll) have made 2D and 
3D hydrodynamic simulations of a slice of SNR, mimicking the 
presence of energetic particles by lowering the adiabatic index 
of the fluid (uniformly in space and time). 

In this work, for the first time we combine all these previ- 
ous approaches; that is, we make full 3D simulations of an SNR 
evolution including a space- and time-dependent model of accel- 
eration and back reaction of particles, to be able to interpret the 
X-ray observations. 
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2. Method 

2.1. Hydrodynamic evolution 

As we are interested in the time-dependent and non-linear in- 
terplay between the SNR and the energetic particles it acceler- 
ates, we resort to numerical simulations . To comp ute the remnant 
evolution, we used the code RAMSES dTevssierll2002h . This 3D 
Cartesian Eulerian code includes a second-order hydrodynamic 
solver, and implements a tree-based structure allowing for ver- 
satile adaptive mesh refinement (AMR). 

Although RAMSES has been already extensively tested and 
used, we had to adap t it to the study of supernova remnants 
(Fraschetti et al. 2009). The main point is that we use a grid that 
is comoving with the contact discontinuity; that is, we work in an 
expanding frame. Because this frame is non-inertial, we have to 
modify the Euler equations. Although it is computationally very 
interesting to factor out the global expansion of the remnant in 
this way, we have to face the numerical instabilities associated to 
quasi-stationary Shockwaves. Still, we can accurately follow the 
SNR evolution as shown in Fig. Q] In the unmodified case, the 
position and spee d of the shocks exactly co incide with analyti- 
cal predictions by Truelove & McKeel(ll999l) . In the simulations 
presented here, we actually simulate one eighth of a sphere, as- 
suming symmetry. 

2.2. Particle acceleration and back reaction 

To compute acceleratio n of particles by s hocks, we used the 
semi- analytical model of Blas](l2002j,IIP04), a non-linear model 
of diffusive shock acceleration (DSA) that takes the back reac- 
tion of accelerated particles on the fluid structure into account. 
This model solves the particle spectrum f(p) and the fluid ve- 
locity profile U(p) jointly as functions of the momentum p of 
particles. It includes the escape of particles with the highest en- 
ergy upstream of the shock, which carry energy away from the 
accelerator. We also include the effect of Alfven wave heating in 
the precursor, which limits shock mo difications, but we d o not 
include magnetic field amplification ( Amato & Blasi 2006). 

As inputs, the model requires (i) information on the shock 
(its speed us and Mach number Ms are provided by the hydro- 
dynamic code, averaged over the remnant surface); (ii) an injec- 
tion recipe (we assume that a fraction r\ of the particles cross- 
ing the shock enter the accelerator, at injection momentum p-^ 
equal to £ times the mean downstream thermal momentum); and 
(iii) a cutoff recipe (we limit the maximum momentum /? max ac- 
cording to the age and the size of the remnant, assuming Bohm 
diffusion, i.e. a diffusion coefficient Dip) = Do p 2 / tJ\ + p 2 
with Do = 3.10 21 cm 2 /s). We consider particle acceleration only 
at the forward shock, as there is less theoretical and observa- 
tional evidence of efficient acceleration at the reverse shock (see 
Ellison et al. 2005 for a discussion of this issue). 

Amongst the outputs, the acceleration model provides the 
total shock compression ratio r lot . To couple the hydrodynamic 
evolution of the remnant with particle a cceleration, we var y the 
adiabatic index of the fluid as done in lEllison et alj d2004l) : at 
each time-step, we compute the index y e ff, which would have 
produced the same ratio r tot , and affect it in RAMSES to the cells 
located just upstream of the shock front. Then y e ff is advected in- 
side the shocked region, so that it remains constant in each fluid 
element, which thus remembers modificati ons induced by parti- 
cle acceleration at the time it was shocked. Elli son et al.l {2004) 
have shown that there is good agreement between such a pseudo- 
fluid approach and two-fluid calculations in ID. 
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Fig. 1. Evolution of the shocks. 

Radial position xs, velocity u$, and Mach number Ms of the forward 
(solid) and reverse (dashed) shocks are plotted versus time, averaged 
over all directions (Mach number is not shown for the reverse shock). 
Blue curves show the unmodified, i.e. purely hydrodynamic case (with 
SNR parameters given in Sect. 12. 3t : red curves show the modified case, 
including the back reaction of accelerated particles (with acceleration 
parameters given in Sect. I3.21 l. 



2.3. Supernova remnant initialisation 

We initialise our simulations a t a young age (h ere 10 years), 
using self-similar profiles from lCheval ier (1983), including the 
pressure of accelerated particles (as computed from the acceler- 
ation model presented in the previous section). Assuming that 
both the ejecta (but for a central uniform core) and the ambi- 
ent medium have power-law density profiles (of indices respec- 
tively n and s), hydrodynamic profiles are obtained by integra- 
tion of ordinary differential equations. 

Here we are interested in a Tycho-like SNR, that is a super- 
nova of type la, referring to 10 51 erg of kinetic energy in 1 .4 solar 
masses, with an n — 7 power-law distribution. We assume a uni- 
form (s = 0) and tenuous (rtH,o =0.1 cm 4 ) ambient medium. 

Finally, we mention that Rayleigh-Taylor instabilities are not 
explicitly seeded in the simulation, but are spontaneously trig- 
gered at the contact discontinuity by numerical perturbations 
seeded by the grid. 
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Fig. 2. Evolution of the remnant structure. 

Top plot: compression factor r of the shock (upper curves show the to- 
tal compression r tot , lower curves show the subshock compression f su b)- 
Middle plot: effective adiabatic index y cff at the shock front. Bottom 
plot: relative positions x of the waves (averaged over all directions). 
From top to bottom, curves correspond to the forward shock, the mixing 
zone boundaries (defined as the region where at least 10% of the density 
is made by a constituent which would not have come here without in- 
stabilities), and the reverse shock. Colour codes the (constant) injection 
fraction rj, rising from blue (almost unmodified case) to red (very modi- 
fied case) asfollows: 10~ 5 , 10~ 4 , 2xl0~ 4 ,4xl0- 4 , 6xl0~ 4 , 8xl0~ 4 , 10~ 3 . 



3. Results 

3.1. Remnant evolution 

The temporal evolution of some key parameters is shown in 
Fig- El The acceleration model depicted in Sect. l2.2l provides the 
compression ratios plotted on top of the figure. For ease of inter- 
pretation, we ran multiple simulations with different fixed injec- 
tion rates 77. If 77 10~ 5 (in dark blue), the system is almost in the 
linear (test-particle) regime, and there is a single strong shock 
of compression ratio r = 4. As we raise 77 to 1CT 3 (as colour 
gets warmer), it progressively enters the non-linear (modified) 
regime, and the shock discontinuity is reduced to a subshock of 
compression r SU b between 3 and 4, whereas the overall compres- 
sion ratio r tot increases to more than 10. (Contrary to ordinary 
hydrodynamics, r tot always significantly depends on the shock 
Mach number, even when the latter is very high, see e.g. Blasi 
120021 ) The corresponding effective adiabatic index y e ff is plotted 
in the middle of the figure. As expected, it decreases from 5/3 




Fig. 3. Maps of the density. 

Snapshot at time 500 years from 3D simulations with formal resolu- 
tion of 1024 3 cells. Left side: slices of the density in the plane z = 0. 
Colour codes phases: ejecta are in green (where their fraction is at least 
10%), ambient medium is in purple (po is its unperturbed mass density). 
Features located just behind the forward shock result from a numeri- 
cal instability. Right side: projection of the square of the ejecta density 
along z-axis. Top half of the figure shows a purely hydrodynamic case, 
bottom half a case includin g the back react ion of accelerated particles 
(with the injection recipe of Blasi et al. 2005 which gives tj ~ 4. 10~ 4 ). 



(thermal fluid) to 4/3 (relativistic fluid) or even below (as par- 
ticles escape). In initially poorly efficient accelerators (low 77), 
shock modifications are small, but steadily increase over time. 
On the other hand, in very efficient accelerators (high 77), back re- 
action effects are very strong even at a very early stage, but then 
decrease over time. In all cases, we see that parameters evolve 
rather slowly after the firs t hundred years. These resu lts are sim- 
ilar to those obtained by iDecourchelle et all J2000) in ID, as- 
suming self-similarity of the hydrodynamic profiles and using a 
different acceleration scheme, which cross-validates the models. 

Finally, the evolution of the relative positions of the waves 
is shown at the bottom of the figure. Two main comments are in 
order. First, the region affected by the Rayleigh-Taylor instabil- 
ity does not seem to be affected by the acceleration of particles. 
Second, the forward shock can get very close to this turbulent 
region if the injection level is high enough (above 5 x 1CT 4 ) 
- the reverse shock also reaches the contact discontinuity, but 
this is caused by the development of the Rayleigh-Taylor in- 
stability, independent of acceler ation efficiency. These o bser- 
va tions agree with the res ults of iBlondin & Ellison! d200ll) and 
of iFraschetti et al.l (120091) . obtained with no space- and time- 
dependent model of acceleration. Our simulations also show that 
the average distance between the forward shock and the con- 
tact discontinuity does not depend much on time in the case of 
efficient acceleration. (Rayleigh-Taylor fingers grow steadily in 
time, but the forward shock moves away as back reaction effects 
are progressively reduced.) 
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3.2. Remnant morphology 

The morphology of the remnant after 500 years is shown in 
Fig. [3] where we compare cases with (top half) and without (bot- 
tom half ) back react i on of particles. We used here the injection 
recipe of Blasi et al.] (120051) . which gives the injection level rj as 
a function of £ (we adopt the common value £ = 3.5) and of the 
subshock compression ratio r su \, (so that injection is switched off 
as the shock gets too modified). Although 77 is time-dependent, 
it remains close to the same value during the whole simulation, 
going down from 4.6 x 10" 4 to 3.9 x 10~ 4 . 

On the left of the figure, we show slices of the density, 
which highlights the remnant's structure. The main effect of par- 
ticle back reaction is apparent, namely that the shocked region 
shrinks. For this level of injection, this causes density to rise 
by a factor typically of two downstream of the forward shock. 
Nevertheless, the size of the region perturbed by the Rayleigh- 
Taylor instability, the development of which depends on the den- 
sity contrast, is basically unchanged. 

On the right of the figure we show projected maps of the 
square of the density of shocked ejecta, which roughly approxi- 
mates the intensity of their thermal X-ray emission. The interior 
of the remnant is filled with a filamen tary texture, structure d over 
scales compatible with the picture of IWarren et al.l (120051) . Even 
in projection, the unstable region is still clearly visible as a bright 
annulus. But here again, cases with and without efficient acceler- 
ation cannot be distinguished from the structure of the shocked 
ejecta alone. 

3.3. Application 

Although the detailed modelling of a specific object is beyond 
the scope of this Letter, we can show the interest of our ap- 
proach by com paring our numerica l results with the observa- 
tional res ults of IWarren et al.l d2005l) regarding Tycho's remnant 
(see also IVolk et al.ll2008l for a study of this object using a ID 
non-linear kinetic model). Their main finding is that the forward 
shock is very close to the contact discontinuity, with a mean 
FS:CD ratio of 1:0.93. Assuming a uniform medium (s = 0, 
a reasonable hypothesis for Tycho), the theoretical ratio is 1 :0.85 
for an ejecta indice n — 7 (assuming a self-similar evolution 
without particle acceleration). As n increases, the predicted ra- 
tio de£reas^s 2 _bu^j^_stillOTly 1^0.89 for n = 14. As pointed out 
by iDwarkadas & Chevalier] dl998l) . an exponential profile may 
be more adapted to the early stages of an SN la; however, this 
would pr oduce an even broader shoc ked region at the time con- 
sidered (Cassam-Chenai et al. 1 12008b . On the other hand, such a 
high FS:CD ratio can be obtained readily with our code through 
the back reac t ion of accelerated particles. The procedure used by 
IWarren et al.l J2005) seems to extract the envelope of the ejecta, 
which can be obtained in our simulations by setting some thresh- 
old on the ejecta fraction - we found the value of 10% used 
here to give comparable results. Then the simulations can repro- 
duce the observ ations provided that 77 is of the order of 5 x 10~ 4 
dVolk et al.l 120081 obtained 77 = 3 x 10~ 4 ). Although we have 
not conducted a complete parametric study yet, we believe that 
uncertainties in SNR parameters or in observed ratios could be 
accounted for by varying the injection level - and reciprocally, a 
good knowledge of the object will allow this crucial parameter 
to be constrained. 

IWarren et al.l d2005l) also found that that the reverse shock 
is deep inside the ejecta, with a mean FS:RS ratio of 1:0.71, 
which is quite puzzling. For n = 7 and s — the predicted ratio 
(without particle acceleration) is 1:0.79. Projection effects might 



cause underestimates of the radius of the reverse shock, but ac- 
cording to our maps this effect is too small to explain the discrep- 
ancy. Considering more realistic profiles from explosion models 
might help understanding the situation. In any case, this obser- 
vation does not favour efficient acceleration at the reverse shock, 
which would shrink the whole remnant's structure even more. 
However, other potentially important effects (like the composi- 
tion and temperature of the ejecta) have to be included in our 
model before drawing firm conclusions. 

4. Conclusion 

We have presented a new code that couples a 3D hydrodynamic 
description of a supernova remnant, allowing consideration of 
hydrodynamic instabilities such as the Rayleigh-Taylor instabil- 
ity, with a realistic kinetic model of non-linear acceleration at the 
blast wave, allowing evaluation of the efficiency of particle ac- 
celeration. We are thus now able for the first time to simulate the 
morphology of SNRs undergoing efficient DSA without limiting 
assumptions on its spatial structure or temporal evolution. 

Our first results confirm the most notable previous findings 
regarding particle back reaction on the SNR morphology: (i) the 
shock structure is all the more compact since acceleration is 
efficient, which provides a clear observational diagnostic; and 
(ii) the development of the Rayleigh-Taylor instability is not sig- 
nificantly affected by acceleration at the forward shock, but it has 
to be taken into account when interpreting observations. 

Regarding the case of Tycho's remnant, comparison of our 
simulations with X-ray observations strengthens the case for ef- 
ficient acceleration of protons at the forward shock. 
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